title: GIS Assignment 2: External Data Sources
packages<-c("cowplot","dismo","leaflet","maps","mapdata","raster","rasterVis","readxl","rgbif","rgdal","tidyverse","utils","rdryad")
sapply(packages, library, character.only=T)
library(ggplot2)
library(maps)
library(OpenStreetMap)
library(mapproj)
library(ggmap)
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
## 
## Attaching package: 'ggmap'
## The following object is masked from 'package:dismo':
## 
##     geocode
## The following object is masked from 'package:cowplot':
## 
##     theme_nothing

Simple distribution map of Chrosomus erythrogaster using gbif data

chrom_dismo <- gbif("chrosomus", species = "erythrogaster", ext = c(-91,-81, 36, 40),
                   geo = TRUE, sp = TRUE, download = TRUE,
                   removeZeros = TRUE)
chrom_dismo_df <- cbind.data.frame(chrom_dismo@coords[,1],
                                  chrom_dismo@coords[,2])
colnames(chrom_dismo_df) <- c("x","y")
us <- map_data("state")
ggplot(data = chrom_dismo_df, aes(x=x, y=y)) +
  geom_polygon(data = us, aes(x=long, y = lat, group = group),
               fill = "white", color="black") +
  geom_point() + xlab("Longitude") + ylab("Latitude") +
  coord_fixed(xlim = c(-93,-76), ylim = c(33, 60)) +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Chromsomus erythrogaster") + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) + 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "lightblue"))

Here I limited my species counts to 2000 individuals per species and added a second species Chrosomus cumberlandensis

chrom_rgbif <- occ_data(scientificName = "Chrosomus erythrogaster",
                 hasCoordinate = TRUE, limit = 2000,
                 decimalLongitude = "-130, -40", 
                 decimalLatitude = "30, 80")

cumb_rgbif <- occ_data(scientificName = "Chrosomus cumberlandensis",
                 hasCoordinate = TRUE, limit = 2000,
                 decimalLongitude = "-130, -40", 
                 decimalLatitude = "30, 80")
Figure 2. Blackside Dace (Chrosomus cumberlandensis)
Figure 3. Southern Redbelly Dace (Chrosomus erythrogaster)
chrom_rgbif_df <- cbind.data.frame(chrom_rgbif$data$species,
                                  chrom_rgbif$data$decimalLatitude,
                                  chrom_rgbif$data$decimalLongitude,
                                  chrom_rgbif$data$stateProvince,
                                  chrom_rgbif$data$year)


cumb_rgbif_df <- cbind.data.frame(cumb_rgbif$data$species,
                                  cumb_rgbif$data$decimalLatitude,
                                  cumb_rgbif$data$decimalLongitude,
                                  cumb_rgbif$data$stateProvince,
                                  cumb_rgbif$data$year)
colnames(chrom_rgbif_df) <- c("species","y","x","state","year")
colnames(cumb_rgbif_df) <- c("species","y","x","state","year")
chrom_rgbif_df <- chrom_rgbif_df[complete.cases(chrom_rgbif_df[1:4]),]
cumb_rgbif_df <- cumb_rgbif_df[complete.cases(cumb_rgbif_df[1:4]),]
ggplot() +
  geom_polygon(data = us, aes(x=long, y = lat, group = group),
               fill = "white", color="black") +
  geom_point(data = chrom_rgbif_df, aes(x=x, y=y, color = species), size = 3) +
  geom_point(data = cumb_rgbif_df, aes(x=x, y=y, color = species), size = 3) +  
  coord_fixed(xlim = c(-107,-65), ylim = c(24,50)) +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Chrosomus erythrogaster and Chrosomus cumberlandensis") + 
  guides(color=guide_legend("Legend", override.aes = list(size = 4))) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) + 
  theme(legend.position = "bottom") +
  theme(legend.title.align = 0.5, legend.box.just = "center") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "lightblue"))

Interactive gbif map that looks at just Chrosomus erythrogaster colored by year recorded.

pal <- colorNumeric(
  palette = "RdYlBu",
  domain = chrom_rgbif_df$year)
colors <- colorFactor(c("#F8766D","#00BA38","#619CFF"), chrom_rgbif_df$year)
leaflet(chrom_rgbif_df) %>% 
  addTiles() %>% 
  addCircleMarkers(chrom_rgbif_df$x,
                   chrom_rgbif_df$y,
                   popup = chrom_rgbif_df$species,
                   label= chrom_rgbif_df$year,
                   color = ~pal(year),
                   weight = 1,
                   fillOpacity = 0.7) %>%
  addMiniMap(position = 'topright',
             width = 100, 
             height = 100,
             toggleDisplay = FALSE) %>%
  addScaleBar(position = "bottomright")%>%
  addLegend("bottomleft", pal = pal, values = chrom_rgbif_df$year,
    title = "Year Recorded",
    labFormat = labelFormat(),
    opacity = 1)
Figure 1. Dung Fly

Dryad data on Dung Flies that looks at clutch size with altitude. Data also collected in Europe, but here we will just look at North America.

Dung_Fly <- read.csv("Dung_Fly.csv")
View(Dung_Fly)
canada <- getData('GADM' , country = "CAN", level=0)
world <- map_data("worldHires")
main_map <- ggplot(Dung_Fly, aes(Longitude, Latitude)) +
  geom_polygon(data = us, aes(x=long, y = lat, group = group),
               fill = "gray", color="white") +
   geom_polygon(data = canada, aes(x=long, y = lat, group = group),
               fill = "gray", color="white") +
  geom_point(aes(color = Altitude, size = `Clutch.Size..First.clutch.`)) +
  coord_fixed(xlim = c(-125,-66), ylim = c(20,60)) +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Dung Fly Clutch Size and Altitude") + 
  guides(color=guide_legend("Altitude", override.aes = list(size = 3))) + 
  guides(size=guide_legend("Clutch.Size..First.clutch.")) +  
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) + 
  theme(legend.position = "right") +
  theme(legend.title.align = 0.5, legend.box.just = "center") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "lightblue"))
main_map

Bioclim data with consideration of a few environmental factors

bioclim <- getData(name = "worldclim", res = 2.5, var = "bio", path = "./")
names(bioclim) <- c("Ann Mean Temp","Mean Diurnal Range","Isothermality","Temperature Seasonality",
                    "Max Temp Warmest Mo","Min Temp Coldest Mo","Ann Temp Range","Mean Temp Wettest Qtr",
                    "Mean Temp Driest Qtr","Mean Temp Warmest Qtr","Mean Temp Coldest Qtr","Annual Precip",
                    "Precip Wettest Mo","Precip Driest Mo","Precip Seasonality","Precip Wettest Qtr",
                    "Precip Driest Qtr","Precip Warmest Qtr","Precip Coldest Qtr")
bio_extent <- extent(x = c(
  min(chrom_rgbif_df$x),
  max(chrom_rgbif_df$x),
  min(chrom_rgbif_df$y),
  max(chrom_rgbif_df$y)))
bioclim_extent <- crop(x = bioclim, y = bio_extent)
bioclim_model <- bioclim(x = bioclim_extent, p = cbind(chrom_rgbif_df$x,chrom_rgbif_df$y))
presence_model <- dismo::predict(object = bioclim_model, 
                                 x = bioclim_extent, 
                                 ext = bio_extent)
gplot(presence_model) + 
  geom_polygon(data = us, aes(x= long, y = lat, group = group),
               fill = "gray", color="black") +
  geom_raster(aes(fill=value)) +
  geom_polygon(data = us, aes(x= long, y = lat, group = group),
               fill = NA, color="black") +
  geom_point(data = chrom_rgbif_df, aes(x = x, y = y), size = 2, color = "black", alpha = 0.5) +
  scale_fill_gradientn(colours=c("brown","yellow","darkgreen"), "Probability") +
  coord_fixed(xlim = c(-107,-65), ylim = c(24,50)) +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Probability of Chrosomus erythrogaster Occurrence") + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) + theme(legend.position = "right") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "lightblue"))

colors <- c("brown","yellow","darkgreen")
leaflet() %>% 
  addTiles() %>%
  addRasterImage(presence_model, colors = colors, opacity = 0.8) %>%
  addCircleMarkers(chrom_rgbif_df$x,
                   chrom_rgbif_df$y,
                   weight = 1,
                   color = "green",
                   fillColor = "blue",
                   fillOpacity = 0.7) %>%
  addMiniMap(position = 'topright',
             width = 100, 
             height = 100,
             toggleDisplay = FALSE) %>%
  addScaleBar(position = "bottomright")
state <- map_data("state")
county <- map_data("county")
point <- data.frame("x" = -84.7471, "y" = 34.6953)
View(chrom_rgbif_df)
inset <- ggplot() + 
  geom_polygon(data = us, aes(x=long, y = lat, group = group),
               fill = "yellow", color="green") +
  geom_polygon(data = canada, aes(x=long, y = lat, group = group),
               fill = "pink", color="green") +
  geom_point(data = Dung_Fly, aes(x=Longitude, y=Latitude, color = Clutch.Size..First.clutch.), size = 5) +
  coord_map(xlim = c(-130, -70), ylim = c(24,50), "polyconic") + 
  theme(panel.background = element_rect(fill = "lightblue"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
        axis.line=element_blank(), axis.text.x=element_blank(), axis.text.y=element_blank(),axis.ticks=element_blank(), 
        axis.title.x=element_blank(), axis.title.y=element_blank()) +
theme(plot.margin=grid::unit(c(0,0,0,0), "mm"))
inset

ggdraw() +
draw_plot(main_map) + 
draw_plot(inset, x = 0.50, y = 0.50, width = 0.50, height = 0.50)